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Abstract 


In this article we consider a new derivation mode for generalized 
communicating P systems (GCPS) corresponding to the functioning of 
population protocols (PP) and based on the sequential derivation mode 
and a fairness condition that permits to ensure a particular sequence of 
configurations. We show that PP can be seen as a particular variant of 
GCPS. We also consider several stochastic evolutions satisfying different 
fairness conditions and particularly focus on those corresponding to the 
run of a Gillespie’s SSA. This permits to further describe the dynamics 
of GCPS by a system of ODEs when the population size goes to the 
infinity. 

Keywords: Population protocols, P systems, Stochastic simulation, 
Gillespie’s algorithm. 


1 Introduction 


The notion of a generalized communicating P system was introduced in 
[25], with the aim of providing a common generalization of various purely 
communicating models in the framework of P systems. 

A generalized communicating P system, or a GCPS for short, corre- 
sponds to a hypergraph where each node is represented by a cell and each 
edge is represented by a rule. Every cell contains a multiset of objects which 
— by communication rules - may move between the cells. The form of a 
communication rule is (a,i)(b,7) > (a,k)(b,1) where a and b are objects 
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and i,j,k,l are labels identifying the input and the output cells. Such a 
rule means that an object a from cell i and an object b from cell 7 move 
synchronously to cell & and cell /, respectively. In this respect, the model 
resembles the Petri Net formalism [20] where tokens from various input 
places come along together to fire a given transition and then fork out to 
destination places, see [25, 5] for more details. 


Depending on the form of the communication rules, several restrictions 
on communication rules can be introduced. Due to the simplicity of their 
rules, the generative power of such restricted systems is of particular interest 
and it has been studied in detail. In [25, 10, 9, 19] it was proved that eight of 
the possible nine restricted variants (modulo symmetry) are able to generate 
any recursively enumerable set of numbers; in the ninth case only finite 
sets of singletons can be obtained. Furthermore, these systems are able to 
achieve this generative power even with relatively small numbers of cells and 
simple underlying (hypergraph) architectures. In [9] a further restriction is 
introduced by considering that the alphabet of objects is a singleton (like in 
Petri Nets) and it is shown that computational completeness can be achieved 
in four of the restricted variants. 


Population protocols (PP) have been introduced in [1] (see [3] for a 
survey) as a model of sensor networks consisting of very limited mobile agents 
with no control over their own movement. A population protocol corresponds 
to acollection of anonymous agents, modeled by finite automata, that interact 
with one another to carry out computations, by updating their states, using 
some rules. Their computational power has been investigated under several 
hypotheses, in most of the cases restricted to finite size populations. In 
particular, predicates stably computable in the original model have been 
characterized as those definable in Presburger arithmetic. The article [6] 
studies the convergence of PP when the population size goes to the infinity. 


The evolution of a PP follows a particular fairness condition: an exe- 
cution is fair if for all configurations C that appear infinitely often in the 
execution, if C is predecessor of a configuration C’, then C’ appears infinitely 
often in the execution. We observe that the fairness condition is used to 
select particular evolutions among all possible runs. In a more general case, 
the fairness was investigated in the area of transition systems and there are 
several possible definitions that can be used to yield the desired behavior [7]. 
We consider such conditions in the case of GCPS systems and obtain a new 
derivation mode which we call fair sequential mode (fs-mode). We further 
study the dynamic behavior of the system in this mode. There are several 
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possible evolution strategies implementing different fairness conditions and 
we are particularly interested in the stochastic ones. We focus on the strategy 
corresponding to a run of the Gillespie stochastic simulation algorithm (SSA). 
Using the correspondence between SSA and ODEs (assuming mass-action 
kinetics) we show that the dynamics of the system can be represented by 
a system of ODEs when the population size goes to the infinity. We also 
consider the converse problem and we give sufficient conditions for a system 
of ODEs to be represented by a GCPS system working in concentration- 
depended stochastic implementation of the fs-mode. We consider several 
examples of GCPS simulating Lotka-Volterra (predator-prey) behavior or 
computing approximations of algebraic numbers. 
We remark that this paper is an extended version of [24]. 


2 Background 


In this section we recall some basic notions and notations used in membrane 
computing, formal language theory and computability theory. For further 
details and information the reader is referred to [18, 19, 21]. 

An alphabet is a finite non-empty set of symbols. For an alphabet V, 
we denote by V* the set of all strings over V, including the empty string, A. 
The length of the string « € V* is the number of symbols which appear in x 
and it is denoted by |x|. The number of occurrences of a symbol a € V in 
x € V* is denoted by |z|,. If x € V* and U C V, then we denote by |z|y 
the number of occurrences of symbols from U in z. 

A finite multiset over V is a mapping M : V —> N; M(a) is said to be 
the multiplicity of a in M (N denotes the set of non-negative integers.) A 
finite multiset / over an alphabet V can be represented by all permutations 
of a string x = aM (a) aN (02) ph Oe Vv", where a, 6 V, 15.9 5 np 
represents M in V*. If no confusion arises, we also may use the customary set 
notation for denoting multisets. The size of a finite multiset M, represented 
by « € V* is defined as Ngey|zla. The difference and the sum of two 
multisets m;, and m2 are denoted as m, — m2 and m, + mz respectively. 


2.1 Population Protocols 


A population protocol is the quintuple P = (Q, %,1,w,6), where Q and © are 
alphabets (elements of Q are called states and those of © input symbols), 
1: 4 + Q is the initial state mapping, w : Q > {0,1} is the individual 
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output function and 6 C Q* is a total transition function. Function 6 can 
be described by using the notation qig2 > qiqh, for (q1,q@2,94,95) € 6. By 
convention we consider that qiq2 > qig2 when no rule involving qq2 in the 
left hand side is present. 

The computation of a population protocol on the input multiset w € * 
is performed by n agents (n = |w|), each of them having a state that is 
an element of Q. Two agents being in states gq, and qo, respectively, may 
interact by the rule qq2 — qq from 6. As a result of such an interaction 
the states of the agents change to qj and qj, respectively. 

More formally, a configuration of a PP is a multiset C’ over Q. A 
transition between two configurations C and C’, written as C > C’, is 
obtained by a single interaction between the agents, i.e., C > C’ if C’ = 
C —{q,q}+{d, 95} for some qi, q2 € C and qigz > gigs € 6. We remark 
that the transition does not change the size of the configuration. 

An execution of the population protocol on the input multiset w is 
defined as an infinite sequence of configurations Co, C 1, C2,... obtained by 
a sequence of transitions C; => Cj41 for all 7 > 0, starting from the initial 
configuration Cp = 1(w). 

An execution is fair if for all configurations C' that appear infinitely 
often in the execution, if C + C’ for some configuration CO’, then C’ appears 
infinitely often in the execution. 

The output value of configuration C, denoted by w(C), is defined as 
follows: 

0, if w(q) = 0, for all gE C, 
gC) =< J. if w(q) = 1, for all gE C, 
undefined, otherwise. 


Let p be a predicate over multisets of elements of . Predicate p can 
be considered as a function whose range is {0,1} and whose domain is the 
collection of these multisets. The predicate p is said to be computed by a 
PP P as defined above if for any input multiset w and every fair execution 
of P on the input w, there exists a number k € N such that w(Cy) = p(w) 
for all k! > k. 

The following result was proved in [1, 2]: 


Theorem 1 A predicate is computable in the population protocol model iff 
it 1s semilinear. 


Recall that semilinear sets are known to correspond to predicates on 
counts of input agents definable in first-order Presburger arithmetic [17]. 
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2.2 Generalized Communicating P Systems 


Next we recall the basic definitions concerning generalized communicating P 
systems [25]. 


Definition 1 A generalized communicating P system (a GCPS) of degree 
n, where n > 1, is an (n+ 4)-tuple II = (O, E,w,...,Wn,R,h) where 


1. O is an alphabet, called the set of objects of II; 
2. E CO; called the set of environmental objects of I; 


38. w,€ OF, 1 <i<n, is the multiset of objects initially associated with 
cell 2; 


4. Ris a finite set of interaction rules (or communication rules) of the 
form (a,i)(b, 7) + (a,k)(b,1), where a,b € O,0 <i,9, k,l <n, and if 
i =0 and j =0, then {a,b} N(O\ E) £0; ie, ag E and/or b ¢ E; 


5. he {1,...,n} is the output cell. 


The system consists of n cells, labeled by natural numbers from 1 to n, 
which contain multisets of objects over O; initially cell 7 contains multiset 
w; (the initial contents of cell 7 is w;). We distinguish an additional special 
cell, labeled by 0, called the environment. The environment contains objects 
of E in an infinite number of copies. 

The cells interact by means of the rules (a,7)(b,7) > (a, k)(b,1), with 
a,b€ O and 0 <i,j,k,1 <n. As the result of the application of the rule, 
object a moves from cell i to cell k and b moves from cell j to cell l. If two 
objects from the environment move to some other cell or cells, then at least 
one of them must not appear in the environment in an infinite number of 
copies. Otherwise, an infinite number of objects can be imported in the 
system in one step. 

A configuration of a GCPS II, as above, is an (n+1)-tuple (zo, 21,.--, Zn) 
with z € (O\ E)* and z € O*, for all 1 <i < n; 2 is the multiset of 
objects present in the environment in a finite number of copies, whereas, for 
alll <i<n, z% is the multiset of objects present inside cell 7. The initial 
configuration of II is the tuple (A, wi,..., Wn). 

Given a multiset of rules R over R and a configuration C = (zo, 21,.--, Zn) 
of II, we say that R is applicable to C if all its elements can be applied si- 
multaneously to the objects of multisets zp, z1,..., Zn, such that every object 
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is used by at most one rule. Then, for a configuration C = (zo, 21,.--, 2n) 

of II, a new configuration C’ = (z@, z{,..-,2/,) is obtained by applying the 

rules of R in a non-deterministic maximally parallel manner: taking an 

applicable multiset of rules # over R such that the application of R results 
/ 


in configuration C’ = (zp, z4,-.-,2,) and there is no other applicable multiset 


of rules R’ over R which properly contains FR. 

It is also possible to replace the maximally parallel strategy of rule 
application by other strategies, called derivation modes (in the context of 
the present paper, the terms mode and strategy are used indifferently). A 
derivation mode lies in the heart of the semantics of P systems and it permits 
to specify which multiset among different possible applicable multisets of 
rules can be applied. When P systems were introduced, only the maximally 
parallel derivation mode was considered which states that corresponding 
multisets should be maximal, i.e., non-extensible. With the appearance of 
the minimal parallel derivation mode [8] the concept of the derivation mode 
had to be precisely defined and [12] presents a framework that permits to 
easily define different derivation modes. 


One application of a multiset of rules satisfying the conditions of a 
derivation mode represents a transition in II from configuration C' to con- 
figuration C’. A transition sequence is said to be a successful generation 
by II if it starts with the initial configuration of II and ends with a halting 
configuration in which a halting condition is satisfied. Generally, the total 
halting condition is considered, being true for a configuration where no 
further transition step can be performed. 


We say that II generates a non-negative integer n if there is a successful 
generation by II such that n is the size of the multiset of objects present 
inside the output cell in the halting configuration. The set of non-negative 
integers generated by a GCPS II in this way is denoted by N(II). It is also 
possible to use GCPS as acceptors, in this case an input multiset is accepted 
if the system halts on it. 


In [25] it is shown that GCPS are able to generate all recursively 
enumerable languages. Moreover this result can be obtained by using various 
restrictions on the type of rules (i.e. induced hypergraph structures), on 
the number of membranes and on the cardinality of the alphabet. We refer 
to [25, 10, 9] for more details. 

If the cardinality of the alphabet O is equal to one, then we refer to 
the corresponding symbol as a token (denoted by e). Hence, we assume that 
O = {e}. We observe that such systems are similar to Petri Nets having a 
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restricted topology. This is especially visible if a graphical notation is used. 
However, the maximal parallelism and the concept of the environment are 
specific to P systems, so we place this study in the latter framework. A 
converse study of P systems from the point of view of Petri Nets can be 
found in [13]. For more details on Petri Nets and membrane computing we 
also refer to [19]. 

In this article we shall consider the dynamics of the configuration of 
GCPS, so we are no more interested in computation (and halting evolutions). 


3 Fair Sequential Derivation Mode 


In this section we are interested in the relation between PP and GCPS. We 
show that in terms of structure PP and GCPS are quite similar, the main 
differences concern the environment and the derivation mode. We define 
a new fair sequential mode (fs-mode) for GCPS and hence we are able to 
encode any PP in a GCPS w.r.t. their dynamics. In some conditions, an 
equivalence between the two models can be stated. 

We consider several fairness conditions for GCPS and their implemen- 
tations with stochastic and Gillespie-like strategies. 


3.1 GCPS in Fair Sequential Mode 


For GCPS, mainly the maximally parallel derivation mode is investigated 
with several attempts to investigate asynchronous or minimally parallel 
derivation mode, see [19] for more details. The derivation mode of PP is 
very particular — it corresponds to a sequential strategy where only one rule 
is applied at each step, like in Petri Nets, but with an additional fairness 
condition corresponding to the definition 3 below. We can consider such a 
strategy in GCPS case as well. 

The notion of fairness was originally developed in the context of transi- 
tion systems to select a particular subset of all the possible executions of 
a system. There exist many ways of defining fairness conditions. We give 
below three different fairness conditions widely found in the literature [7]. 


Definition 2 A sequence of transitions Co > C, => ... 1s strong locally fair 
(SLE) iff for every rule r, if r is applicable in infinitely many configurations 
C;, then there are infinitely many j such that one passes from Cj to Cj+1 
(Cj => Cj41) by applying rule r. 
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Definition 3 A sequence of transitions Cop > Ci >... is strong globally 
fair (SGF) iff for every C and C" such that C > C’, if C = C; for infinitely 
many i, then C; = C" for infinitely many j. 


Definition 4 A sequence of transitions Co > C, =>... is weakly fair (WF) 
iff for every rule r, if it exists 7 such that r is applicable in configuration 
C; for alli > 7, then for infinitely many k one passes from Cy to Cy41 
(Cy => Cyii) by applying rule r. 


When assuming that the number of configurations is finite (this is the 
case for classical PP for example), the definition 3 can be easily rephrased 
as follows: a computation Cp > C, =>... is fair if 


e there exists a non-negative integer N such that configuration Cy 
belongs to a terminal strongly connected component of the state 
graph!; and 


e any state of this terminal strongly connected component appears 
infinitely often in the execution. 


In the context of this paper, we are interested in the study of the 
behavior of GCPS in sequential mode with some fairness condition: we 
say that a GCPS II evolve in a fair sequential derivation mode (fs-mode) 
w.r.t. to the fairness condition F if the sequence of transitions between 
configurations of IT fulfills the fairness condition F. 


Encoding PP in GCPS. It can be easily seen that both PP and GCPS 
are particular instances of multiset rewriting. Indeed, in both cases the 
underlying data structure is multiset (obtained in a direct way for PP and 
by attaching the indices of membranes to the objects in GCPS) and the 
evolution rules clearly correspond to multiset rewriting rules with both left 
hand and right hand sides of size two. So, the translation of a PP to a one- 
symbol GCPS can be easily done as follows. Given a PP P with set of states 
Q (for convenience we suppose that Q = {1,...,n}) and transition relation 6 
in an initial configuration Co, let define GCPS IIp = (O, EF, w1,..., Wn, R, 1) 
as follows 


© 0=E={e}, 


‘In this directed graph, nodes correspond to the configurations C, and two nodes C 
and C” are directly linked if C > C’. 
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e w= ee k= |Colq for any q € Q, 


e R={(e,q1)(©,q2) — (©, 4) (©, 95) | 192 > Ody € OF. 


The above system encodes each state g of P by a token e present in membrane 
labeled by q. It is easy to see that there is a one-to-one correspondence 
between the sequences of transitions of P and of Ip. 

Conversely, given a GCPS II with an empty environmental set (i.e., EH = 
() and with a deterministic rule set R (eventually extended to a total function 
with the identity), it is possible to define PP Py with Q = O x {1,...,n} 
and 6 = R, 1.e., each couple (object, cell) corresponds to a state in Py. The 
initial configuration Co € Q* of Pr is given such that |Co|(o,:) = |wilo, for 
all (0,1) € Q. 

Such an encoding works only if the set EF is empty, otherwise it is not 
possible to represent an infinite multiset (corresponding to the environment) 
using PP. 

From these considerations, it is trivial to establish that PP = GCPS 
in fs-mode with respect to SGF when E = 0 and R is deterministic. 

If we consider an encoding function 1 as defined for PP and the halting 
condition corresponding to the stabilization of the w-image of the configura- 
tion, then as an immediate consequence of [1, 2], we obtain that any GCPS 
working in fs-mode with respect to SGF and that has a deterministic rule 
set not involving the environment can only accept semilinear sets. 

Conversely, we also obtain that any PP working in maximally parallel 
mode (i.e., a maximally parallel number of interactions can happen at each 
step) is computationally complete if the number agents in some particular 
state go is unbounded (going to the infinity). 


3.2 Simulation of GCPS in FS-Mode 


When simulating a transition system the notion of fairness should be replaced 
by a local strategy that guarantees the validity of the condition. 

Simulating a GCPS in fs-mode consists in defining an algorithm, called a 
scheduler, that decides at each time step which applicable rule of the current 
state has to be applied. A scheduler can be seen as a local strategy that 
guarantees the validity of the fairness condition during the whole execution. 
It is characterized by the subset of sequences of transitions it can generate 
and then by the fairness condition it respects. 

Among all possible implementations of schedulers, Markovian processes 
feature prominently since they do not require any history or global knowledge 
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on the state space, and they provide a modeling tool useful in many domains 
(like in the simulation of population behaviors or in distributed algorithmics). 
Such a Markovian process corresponds to a labeling of the state graph arrows 
C =?r(©) C’ by astatic probability p,(C) that only depends on configuration 
C and on the applied rule r = (¢, qi) (¢, q2) > (¢, ¢))(¢,q4) (without loss of 
generality we can focus on one token GCPS). We give here two examples of 
Markovian processes: 


1. Equiprobable implementation: C >1/'c C’ where kc denotes the car- 
dinality of the set {C’|C > C’}. 


2. Concentration-dependent implementation: C =?r©) C’ where pr(C) is 
proportional to h,(C), the number of distinct combinations of tokens 
in C that activate r and to a stochastic constant c, that only depends 
on r. The number h,.(C) is given by 


IClalCla ifm £0,q2 40,91 4 Ge 
_ J IClad€la -1) ifa 40,0 = 4 
ne) 7 Cla if qa =0 (1) 
IC las if a= 0 


The two last cases hold when the environment (containing an infinite 
number of tokens) is involved in the rule. 


It is possible to give sufficient conditions for a Markov process to 
generate fair sequences of transitions with probability 1 with respect to some 
fairness condition. For example, Theorem (1) from [4] establishes that any 
bounded fair transition system, that is a transition system where there exists 
a real number ¢ with 0 < ¢ < 1 such that p, > c for all transitions C >? C’, 
generates fair sequences of transitions with probability 1 w.r.t. SLF. This 
result can be straightforwardly applied for the equiprobable implementation: 
since any GCPS has sets of rules R of finite size, we have 1/kc > 1/|R| for 
any configuration C’. Therefore, the equiprobable implementation generates 
fair sequences of transitions with probability 1 w.r.t. SLF. 

To our knowledge, there is no such general result applicable for the 
concentration-dependent implementation. Nevertheless, for the class of 
GCPS equivalent to PP (see Section 3.1), a random pairing scheduler (con- 
sisting in choosing an ordered pair of agents at random, independently and 
uniformly from all pairs) generates fair sequences of transitions with proba- 
bility 1 for SGF [1]. The proof of this result relies on the fact that a PP has a 
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finite state graph and that the associated Markov chain leaves non-terminal 
strongly connected components and visits infinitely often all states of one of 
the terminal strongly connected components with probability 1. Following 
the same reasoning, this result also holds for any concentration-dependent 
scheduler on GCPS with an empty environment. In fact, the usual PP 
random pairing scheduler is a particular concentration-dependent scheduler 


where p,(C) = reo: 


3.3. Gillespie’s SSA 


A usual abstraction in the simulation of biochemical systems consists in 
considering the system (e.g., a bacterium) as a homogeneous chemical solution 
where the reactions of the model are taking place. D.T. Gillespie has proposed 
in [14] a stochastic simulation algorithm (SSA) for producing the trajectories 
of such a chemical system by computing the nezt reaction and the elapsed 
time since last reaction occurred. In the following, we show that Gillpesie’s 
SSA provides a concentration-dependent scheduler where the transition 
probabilities depend on a stochastic coefficient associated with each rule. 


Specification of the SSA. Let yz be a chemical reaction. The probability 
that ys takes place during an infinitesimal time step is proportional to: 
e c,, the stochastic reaction constant? of reaction 1; 


e h,(S), the number of distinct molecular combinations that can activate 
reaction jy; it depends on the current chemical state S; 


e dr, the length of the time interval. 


Gillespie proved that the probability P(7, |S)dr that, being in a chemical 
state S, the next reaction will be of type yu and will occur in the time interval 
(¢+7,t+7+dr) is: 


P(r, u|S)dr = a,(S) e U5) T dr 


The function a,,(S) = c,hy(S) is called the propensity of reaction j, and 
a(S) = >>), a,(S) is the combined propensity of all reactions. 

This probability leads to the first straightforward Gillespie’s exact 
stochastic simulation algorithm called the first reaction method. From a 


?Evaluating the stochastic constants is one of the key issues in stochastic simulations 
of biochemical reactions. 
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chemical state S, it consists in choosing an elapsed time 7 for each reaction py 
according to the probability P(7, |S). The reaction with the lowest elapsed 
time is selected and applied on the system making its state evolve. A new 
probability distribution is then computed for this new state and the process 
is iterated. 


Gillespie’s Simulation of a GCPS. Géillespie’s SSA can be used in 
wider range of contexts than only (bio)chemical modeling. Applied to a 
GCPS II, SSA gives a way to simulate a continuous-time Markov chain with 
the states corresponding to the configurations of II and with transitions 
between states corresponding to a single application of a rule of II. From this 
point of view, the definition of a GCPS is extended to associate with each 
rule r : (a,7)(b, 7) 3° (a,1)(b,k) a stochastic coefficient c, (t.e., a positive 
real number) that characterizes the kinetics of the rule in the Markov 
chain. The probability for a transition between two configurations C’ and C’ 
corresponding to the application of a rule r with stochastic coefficient c, is 
given by 
Ma ar(C) = Crhy(C) 
a(C) a(C) 

that directly corresponds to the definition of a concentration-dependent 
scheduler. 

As an immediate consequence, Gillespie’s algorithm generates fair ex- 
ecutions of PP w.r.t. SGF; in other words, it allows the computation of 
semilinear predicates with probability 1. 

The use of Gillespie’s SSA allows us to deal with interesting population 
behaviors in GCPS. It is illustrated in the two next sections where the 
Gillespie’s simulation of GCPS is used to characterized population dynamics. 


4 Application to the Lotka-Volterra Model 


A paradigmatic example illustrating how GCPS allows a well suited specifica- 
tion of population behaviors consists of the description of a process inspired 
by the Lotka-Volterra model. 


4.1 The Lotka-Volterra Model 


The Lotka-Volterra process was introduced by Lotka as a model of coupled 
auto-catalytic chemical reactions, and was investigated by Volterra as a 
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model for studying an ecosystem of predators and preys [11]. This model 
specifies how two coupled populations (of chemicals or individuals) Y; (the 
preys) and Y2 (the predators) behave. It may be presented using the following 


ODEs 
dY, dY> 


dt dt 
In [14], D.T. Gillespie proposes the study of this model from a discrete and 
stochastic point view based on the system of chemical reactions 


= (cy = c2Y2) Yi = (c2Y1 = c3)Y2 (2) 


Cl co C3 
Y, —— 2Y; Y, + Yo ———> 2Yo Y, ———.. (3) 


The dynamics of these reactions is conveniently characterized using the 
predator-prey interpretation. The first rule states that a prey Y; reproduces. 
The second rule states that a predator Y2 reproduces after feeding on prey 
Y,. Finally, the last rule specifies that predators Y2 die of natural causes. 
Coefficients c; are the rates of the three reactions. The correspondence 
between the two models relies on the fact that the trajectories of Gillespie’s 
SSA tend to the solutions of the ODEs system given by the law of mass 
action on the reactions when the number of chemicals tends to the infinity 
(see Section 5). This result is due to the particular application of Kurtz’s 
theorem [15] to chemical systems. 


4.2 Lotka-Volterra GCPS Definition 


The model above does not fulfill the GCPS requirements since the first 
and last reactions are not pairwise interactions. We propose to extend 
reactions (3) by considering a renewable resource X for Y; as a third species?: 
the molecular level of X remains constant whatever its production or its 
consumption. The extended system of reactions is: 


J / / 
= Cy Co _ & _ 
Kay —+.27y Ye y—_3 2%, Yotx—— ox 4) 


The use of a pairwise interaction in the last reaction can be interpreted as 
a competition between the two predator behaviors: a predator in presence 
of preys eats and reproduces (second reaction); a predator in absence of 
prey (represented by the grass) dies (third reaction). Moreover, with the 


3We use the same notation as in [14] to express that the food resource X is assumed 
renewable. 
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Figure 1: Runs of the Lotka-Volterra model with renewable (left column) and non- 
renewable (right column) resources for different initial states (kinetics rates equal 1): 
Y, = 500, Y2 = 1500, X = X = 1000 (first row), Yi = Yo = X = X = 1000 (second row), 
Y; = Y2 = 1000, X = X = 10000 (third row), Y: = 1000, Y2 = 0, X = X = 10000 (fourth 
row). The solid red line represents preys, the dashed green line predators, and the dotted 
blue line resources. The two first rows show that both dynamics exhibit the same properties 
as presented in [14] (particularly, in second row, oscillations araise from an equilibrium 
initial state for the ODEs). The third row shows the difference in the dynamics when the 
resource size is ten times larger than the population size. The last row shows the difference 
in the dynamics when the predator population is empty. The simulations have been done 
using the general simulation language MGS (http://mgs.spatial-computing.org) that 
allows an easy implementation of all models of the present article [22, 16, 23]. 
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hypothesis that the number of X remains constant, systems (3) and (4) 
have equivalent dynamics in terms of propensity when c; = Xc}, c2 = & 
and cz; = Xc3: propensity of the first reaction of system (3) is given by 
oe Oy ee ae Sime cot X Y;, that is the expected propensity of the first 
reaction of system (4); second reactions of both systems are the same; a 
similar reasoning can be done for the last case. 


System (4) is only composed of pairwise interactions that satisfy con- 
dition 4 of Definition 1. Thus, it can be easily expressed as a one-symbol 
GCPS, denoted Izy, with rules R: 


(©, 0)(@,1) +% (¢,1)(¢,1) 
)(,2) +% (0, 2)(@, 2) 
0 


—_~ > 
e ..°e 
NO FR 


)(e,0) +% (¢,0)(¢, 0) 


where membrane indices 0, 1 and 2 represent the environment (an infinite 
source of X), the preys Y; and the predators Y2, respectively. 


Let now consider the previously defined concentration-dependent sched- 
uler. A transition C >? C’ has a probability of p, = a,/ao with the 
propensity function a, = c,h: ¢, is the rate of the corresponding reaction 
in (3) and h, is given by equation (1) accordingly to Izy. The reader is 
invited to pay attention that even if the environment is an infinite source 
of X (instead of a constant one), the dynamics is well taken into account: 
rules involving environmental objects have probabilities that do not depend 
on the number of these objects, see equation (1). In this respect, any com- 
putation of Izy represents a run of the Gillespie’s SSA of reactions (4). As 
a consequence, Izy is an exact model of the original Lotka-Volterra system. 


It has to be remarked that Izy cannot be described by any PP since the 
environment objects are involved in its definition. A possible specification of 
the Lotka-Volterra equations may be obtained within a PP by considering 
X as a non-renewable resource. Such a definition has been realized (taking 
reactions (4) and substituting X by X.) However, due to the limitation of 
resource, this system does not respect the original Lotka-Volterra system dy- 
namics anymore. For example, without any predators, a population of preys 
stabilizes in this model, while in the original model it grows exponentially. 
Figure 1 gives some examples of simulations of the Lotka-Volterra model 
considering renewable and non-renewable resources. 
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5 GCPS With Infinite Population Size 


In this section, we study how the dynamics of GCPS with infinite population 
can be described by a system of ODEs. 


5.1 GCPS Population Dynamics 


As previously said in the presentation of the Lotka-Volterra model, the mass 
action law gives a way to associate a system of ODEs with a system of 
chemical reactions. It is possible to reverse this method and to give a GCPS 
system whose population dynamics (when the population size tends to the 
infinity) will correspond to some dynamics given by a system of ODEs. 
Let us consider the ODEs system defined on set of variables {Y1,..., Ya } 


of the form ie 
ie = SC a Yi¥n — Do (big + by) VY) (5) 
ik J 


where coefficients ay and b;; satisfy the following conditions: 
1. for all i, j,k, ap > 0 and b;; > 0; 


2. for all j,k such that bj, 4 0, there exists either one index % such that 
ais, = 2b;,, or two distinct indices i; and ig such that ap = ai, =U, 
for any other index i, 7 4 ig or (i £ 7, and 7 F¥ ig), ai, = 0. 
The above conditions are sufficient to ensure that >°; aM = 0. Then there 
exists a concentration-dependent implementation of fs-mode GCPS with 
respect to the fairness condition SGF and without rules involving the envi- 
ronment (7.e., a PP) whose behavior is exactly described by ODEs (5) when 
the population size goes to the infinity. Indeed, these equations correspond 
to the mass-action law of a set of rules such that for any j,k with bj, 4 0 


(0, j)(¢,k) Pik (#, i0)(@, zo) or (e, j)(¢, k) —ybik (@, 11)(¢, i2) 


according to the considered possibility of the above condition 2. The reader 
is invited to pay attention to the fact that these equations correspond to 
a wider range of dynamics than the dynamics of second-order chemical 
reactions with two products since they allow the specification of ordered 
interactions (e.g., involving a sender and a receiver as considered in the PP 
literature). This property also holds in PP and suggests that equations (5) 


Generalized Communicating P Systems in Fair Sequential Mode 243 


exactly describe PP dynamics when the size of the populations tends to the 
infinity. 

It is obvious that a more general class of population behaviors is captured 
by concentration-dependent fs-mode GCPS since they do not need to be 
conservative thanks to the environment. Following the idea of equivalence 
between systems (3) and (4) in terms of dynamics, equations (5) can be 
extended with the introduction of a renewable variable Yo. 


5.2 Computational Properties 


We focus here on the original use of PP as a computational model of 
algebraic numbers proposed in [6]. This article investigates the case where 
the computation is independent of the initial contents of the system. 

Using the GCPS terminology, the main idea of [6] is to consider the 
result of a computation as a ratio between the number of tokens in certain 
membrane and the total number of tokens (without taking care of the 
environment) when the population size goes to the infinity and when the 
state of the system converges. The proposed work relies on the definition of 
a particular strategy of execution of the PP: a step of execution consists in 
sampling uniformly and independently of the past two distinct tokens in the 
membrane and let them interact in a sequential mode. The authors of the 
aforementioned article studied the Markov chain associated with PP and 
proved its equivalence to some system of ODEs at the limit. 

We remark that the same kind of result directly arises from considera- 
tions of Section 5.1 since this computational model is captured by one-symbol 
GCPS working in fs-mode with Gillespie concentration-dependent imple- 
mentation. Indeed, the above execution strategy exactly corresponds to a 
Gillespie’s SSA run where the stochastic constants equal 1 for all rules. Thus, 
the study of the model corresponds to the investigation of the sensibility of 
the associated ODEs system. Let us illustrate this point by considering the 
running example of [6] 


(¢, p)(e,p) — (¢,p)(*,m) 
(©, p)(e,m) — (@,p)(©,p) 
(e,m)(e,p) + (©, p)(*,p) 
(e, m)(e,m) — (e,p)(e,m) 


where symbols p and m identify two membranes. It has been shown that 


the ratio where p (resp. m) is the size of the membrane p (resp. m), 


Pp 
p+m? 
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converges to a when the population size goes to infinity. Accordingly to 


equations (5), we associate ODEs with this GCPS as follows 

Hw Xm + 2¥p¥in — ¥p Hw Xn — 2¥p¥m + ¥p 
The stable states of this system are obtained when the two equations vanish, 
that is, when either Yn = —(/2+1)Yp or Yn = (V2 —1)Yp. The first 
solution is incoherent since it involves a negative size of population. The 
second solution trivially leads to the expected result aon = won 

Using this method, we found the 17 algebraic numbers computable 
with the 81 PP involving only two states: {0,2 — V3, a5, 7 pV aT. V2 


1 V17-3 1 V5-1 1 1 2 
aa 17RD ah vai Wee V2, 3) 3g V3 1,1}. In the general 


case, there are ° PP with n states and we do not know how many 
algebraic numbers are computable by them. We conjecture that only solutions 
of the system of equations an = 0 derived from the ODEs system 5 can be 
computed. 


C2.) n 


6 Conclusions 


In this article we investigated connections between population protocols 
and generalized communicating P systems. The two models share the same 
multiset structure and the same type of rules. Traditionally PP are used to 
study population dynamics in the context of distributed algorithmics while 
GCPS are investigated for computational properties. 

By incorporating the derivation mode from PP into GCPS framework 
we obtained a strict inclusion of PP in GCPS working in fs-mode. We then 
considered different fairness conditions and we concentrated on a particular 
implementation of the SGF corresponding to a run of Gillespie’s SSA. As 
a result we obtained that the dynamics of such systems can be described 
by the corresponding system of differential equations. Different questions 
could then be explored, like the investigation of the conditions ensuring that 
the system reaches a stable state regardless of its initial state or ensuring 
that a stable state is never reached for any initial configuration. GCPS 
are in this sense easier to handle than PP because of the environment that 
permits to easily simulate the equivalent of creation or degradation reactions. 
Section 5.1 also considers the converse problem of the construction of a 
GCPS system exhibiting a particular behavior given by a systems of ODEs. 
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It would be interesting to see if the given sufficient conditions are also 
necessary. A mathematical challenge resulting from Section 5.2 is whether 
for any algebraic number « € [0..1] being a root of a system of polynomials 
of degree 2 there is a concentration-dependent evolution implementation of 
GCPS working in the fs-mode with respect to SGF that converges to x. 

We remark that the presented results hold only in the concentration- 
dependent scheduler. By taking an equiprobable scheduler the results are 
different. 

Since Petri Nets can be seen as multiset rewriting, it is clear that the 
results of this paper can be translated to this domain (for Petri Nets with 
specific types of rules and an additional fairness strategy). 

We think that the fs-mode has interesting properties that should be 
further explored. The fairness condition can express a local stochastic 
evolution, so it could be preferable to consider this condition instead of a 
stochastic behavior. Another interesting property of the proposed stochastic 
implementation is that Gillespie’s SSA introduces an explicit continuous 
time and discrete events in the model, which do not appear in a GCPS 
description. 
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